Ab-initio theory of superconductivity — II: Application to elemental metals 



M. A. L. Marques, M. Luders,^'^ N.N. Lathiotakis/-^ G. Profeta," 
A. Floris/'^ L. Fast,6-2 A. Continenza,'^ E. K. U. Gross,!^^ ^nd S. Massidda^' * 

'institut fur Theoretische Physik, Freie Universitat Berlin, Arnimallee I4, D-14-195 Berlin, Germany 
^Institut fiir Theoretische Physik, Universitat Wiirzburg, Am Hubland, D-97074 Wiirzburg, Germany 
' Daresbury Laboratory, Warrington WA4 iAD, United Kingdom 
^CASTI - Istituto Nazionale Fisica della Materia (INFM) and Dipartimento di Fisica, 
Universita degli studi dell'Aquila, 1-67010 Coppito (L'Aquila) Italy 
^INFM SLAGS, Sardinian Laboratory for Gomputational Materials Science and Dipartimento di Scienze Fisiche, 
Universita degli Studi di Gagliari, S.P. Monserrato-Sestu km 0.700, 1-09124 Monserrato (Cagliari), Italy 
^ SP Swedish National Testing and Research Institute, P.O.B. 857, S-501 15 Boras, Sweden 

The density functional theory for superconductors developed in the preceding article is applied to 
the calculation of superconducting properties of several elemental metals. In particular, we present 
results for the transition temperature, for the gap at zero temperature, and for thermodynamic 
properties like the specific heat. We obtain an unprecedented agreement with experimental results. 
Superconductors both with strong and weak electron-phonon coupling are equally well described. 
This demonstrates that, as far as conventional superconductivity is concerned, the first-principles 
prediction of superconducting properties is feasible. 

PACS numbers: 74.25.Jb, 74.25.Kc, 74.20.-z, 74.70.Ad, 71.15.Mb 



I. INTRODUCTION 

The recent discovery of superconductivity at around 
40 K in MgB2^ has renewed the attention of the scien- 
tific community on this field. MgB2 is just one in a 
long list of materials that are found to be superconduct- 
ing. This list includes several elemental metals, heavy 
fermion compounds, high-Tc ceramics^, fuUerenes doped 
with alkali atoms'^, etc. The mechanism responsible for 
superconductivity can have different origins. For exam- 
ple, in the elemental metals, in the fuUerenes'^, and also in 
MgB2''^'^ the electron-phonon interaction is the responsi- 
ble for the binding of the Cooper pairs. This situation is 
usually referred to as "conventional superconductivity" . 
On the other hand, it is generally believed that the very 
high transition temperatures exhibited by the high-Tc 's 
are (at least partly) due to Coulombic effects. Following 
the remarkable experimental discoveries of the last years, 
there were numerous theoretical developments, that have 
greatly improved our description and understanding of 
superconductivity. However, the prediction of material- 
specific properties of superconductors still remains one of 
the great challenges of modern condensed-matter theory. 

In this work, we present an ab-initio theory to de- 
scribe the superconducting state. It is based on a density 
functional formulation, and is capable of describing both 
weakly and strongly coupled superconductors. This is 
achieved by treating the electron-phonon and Coulomb 
interactions on the same footing. The main equation 
of this theory resembles the gap equation of the theory 
of Bardeen, Cooper and Schrieffer^. It is, however, free 
of any adjustable parameter and contains effects origi- 
nating from the retarded nature of the electron-phonon 
interaction. The theoretical foundations of our approach 
are presented in the preceding paper^, which will hence- 
forth be referred to as I. As in ordinary density functional 



theory (DFT)^"^°, the complexities of the many-body 
problem are included in an exchange-correlation func- 
tional. In I we use Kohn-Sham perturbation theory^^ 
to derive several approximations for this quantity. In 
the present paper we describe the implementation of our 
theory, and its application to the calculation of supercon- 
ducting properties of elemental metals. The systems un- 
der consideration range from weak-coupling (Mo, Al, Ta) 
to strong-coupling (Nb, Pb) superconductors. By study- 
ing these well-known systems, we illustrate the usefulness 
and accuracy of DFT for superconductors. Furthermore, 
our results serve as a justification for the choices and 
approximations made in I. Further applications of our 
approach to more complex systems like MgB2^ or solids 
under pressure will be presented in separate publications. 

This paper is organized as follows: In Sect. II we give 
a brief summary of the theoretical foundations of our 
work. The exchange-correlation functionals are described 
in Section III. We present three different levels of ap- 
proximation: functionals that retain the full dependence 
on the wave- vector, energy averaged functionals, and hy- 
brid schemes. We proceed with an account of the com- 
putational details of our numerical implementation. In 
Sect. V we present and discuss numerical results obtained 
for the elemental metals. These results include calcula- 
tions of the transition temperature Tc, the gap at zero 
temperature Aq, and thermodynamic properties like the 
specific heat. The last section is devoted to the conclu- 
sions. 



II. METHOD 

In this section we give a brief account of the theoretical 
foundations of DFT for the superconducting state. For 
an in depth description of this theory, we refer the reader 
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to I. 

A correct description of the superconducting state 
has to include the effects of the electron-electron and 
electron-phonon interactions. DFT for superconductors 
is a theory designed to treat on the same footing both 
electronic correlations and the electron-phonon coupling. 
To achieve this unified description, we start with the 
full electron-nuclear Hamiltonian. A multi-component 
is then established using a set of three densities: 
i) the normal electronic density n{r); ii) the anomalous 
density xC^j^"')) which is the order parameter of the su- 
perconducting state; and iii) the diagonal of the nuclear 
density matrix r(^), where ^ is a shorthand for the 
A'' nuclear coordinates {Ri, R2, ■ ■ ■ , Rn}- An extension 
of the Hohenberg-Kohn theorem^'^'^ guarantees a one- 
to-one correspondence between the set of the densities 
{'^C^)) x(^) r(^)} in thermal equilibrium and the set 
of their conjugate potentials. As a consequence, all ob- 
servables are functional of this set of densities. 

We then construct a Kohn-Sham system^ composed of 
non-interacting (but superconducting) electrons and nu- 
clei. The latter interact with each other through an A''- 
body potential, but do not interact with the electrons. 
The Kohn-Sham system is chosen such that the Kohn- 
Sham densities are equal to the densities of the interact- 
ing system. The Euler-Lagrange equations for the Kohn- 
Sham system lead to a set of three coupled equations, 
one of which describing the nuclear degrees of freedom, 
and the other two describing the electrons. These three 
equations have to be solved self-consistently. The nuclear 
equation describes a set of N nuclei under the influence of 
an effective A^-body potential Vg[n, x, ^]{R)^ and has the 
same structure as the usual Born-Oppenheimer equation 
for the nuclei. In this article we are interested in solids 
at relatively low temperature, where the nuclei perform 
small oscillations around their equilibrium positions. In 
this case, we can expand v^[n,x,^]{R) in a Taylor se- 
ries around the equilibrium positions, and transform the 
nuclear degrees of freedom into collective (phonon) coor- 
dinates. The Kohn-Sham electrons obey a system of two 
coupled equations 
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Unk{r) + j d?r' A^{r,r')vnk{r') 

= EnkU„k{r) (la) 

Vnk{r) + jd?r' ^t{ry)unk{r') 

= EnkVnkir). (lb) 



Note that for a periodic solid, we can classify the cor- 
responding solutions according to the symmetry of the 
translational group of the crystal. Thus we label the 
eigenstates with a principal quantum number n and a 
wave-vector quantum number k. The Eqs. (1) have the 
same structure as the Bogoliubov-de Gennes equations, 
and describe a system of non-interacting, superconduct- 
ing, electrons moving under the influence of the effective 
potential v°{r) and the effective pairing field As{r,r'). 



As it is usual in DFT, the effective potentials {v°, Ag, w"} 
are functionals of the set of densities {n, x, F} and include 
Hartree and exchange-correlation contributions. These 
latter terms, the exchange-correlation potentials, are de- 
fined as functional derivatives with respect to the densi- 
ties of the exchange-correlation contribution to the free 
energy F^c- 

The full self-consistent solution of the three Kohn- 
Sham equations is a highly demanding task. For the 
time being, we resort to some additional approximations 
in order to simplify the problem. First, we neglect the de- 
pendence of the nuclear potential on x, which amounts 
to neglecting the effect of the superconducting pair po- 
tential on the phonon dispersion. This effect has been 
measured experimentally^'*, and it turns out to be quite 
small. Furthermore, we assume that the nuclear poten- 
tial w" is well approximated by the Born-Oppenheimer 
potential. Within these approximations, the phonon fre- 
quencies and the electron-phonon coupling constants can 
be calculated using standard density functional linear- 
response methods*^. 

A direct solution of the Kohn-Sham Bogoliubov-de 
Gennes equations (1)*^ is faced with the problem that 
one needs extremely high accuracy to resolve the super- 
conducting energy scale, which is about three orders of 
magnitude smaller than typical electronic energies. At 
the same time, one has to cover the whole energy range 
of the electronic band structure. The problem can be sim- 
plified through the decoupling approximation*^. First we 
assume that w|(r) does not depend significantly on the 
anomalous density x- In fact, we expect that correc- 
tions to Vg{r) are of the order of |xP and therefore much 
smaller than the typical electronic energy scale. The 
normal-state Kohn-Sham eigenfunctions and eigenener- 
gies can then be used to construct approximations for 
the eigenstates of the superconducting phase 

Unk{r) ^ Unk^nkir) ; Vnkir) pa Vnk^Pnkir) ■ (2) 

The ifnk^s are the solutions of the normal-state Kohn- 
Sham equations for band n and wave vector k and can 
thus be calculated using standard electronic structure 
methods. If A/'b is the number of bands, the decoupling 
approximation transforms, for every k point in the Bril- 
louin zone, the 2A^b x 2Ab eigenvalue problem given by 
Eq. (1) into A^b(2 x 2) secular equations. The Kohn-Sham 
eigenenergies then become 



Enk = Jenk + \^ 
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where ^nk are the normal-state Kohn-Sham eigenvalues 
measured relative to the chemical potential ^nk = enk—^- 
The matrix elements A„fc, which are a central quantity 
in our formalism, are defined as 

Ank = j A\j d?r' ^lk{r)A,{ry)^nk{r') . (4) 

Within the decoupling approximation, A„fc is deter- 
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mined by the integral equation 



A„fe = — ^„,feA 



2 ^ ^ ^nk,n'k' 

n' k' 



tanh 



En'k' 



(5) 

where /3 is the inverse temperature. The quantities 2nk 
and K,nk,n'k' are functionals of the pair potential A„fc 
and of the chemical potential /i. Wc note that, although 
the gap equation (5) is static (i.e. it does not depend ex- 
plicitly on the frequency), it includes retardation effects 
through the functionals Znk and ICnk,n'k'- We will come 
back to this point later in the discussion of our results. 

It is possible to view the decoupling approximation 
from a different perspective: It can be shown that the ma- 
trix elements of the pair potential As(r, r') arc diagonal 
with respect to all symmetry-related quantum numbers, 
in particular the Bloch wave-vector k, but also labels 
of the point-group symmetry (which usually are not ex- 
plicitly given). The decoupling approximation amounts 
to neglecting the matrix elements which are off-diagonal 
with respect to the band index n. For a given fe-point in 
the Brillouin-zone, the corresponding states are in gen- 
eral energetically far apart, which justifies the neglect of 
these matrix elements. The only situation where this ap- 
proximation may break down is when two bands of the 
same symmetry cross in the vicinity of the Fermi surface. 
We note in passing that the decoupling approximation is 
(trivially) exact for the uniform superconducting electron 
gas. 

The decoupling approximation can further be justified 
by the fact that a perturbative treatment of the neglected 
off-diagonal terms does not give any contribution in first 
order. To see this, we consider the neglected part of the 
pair-potential 

A(r,r')= Yl 'Pnk{r)Ak,nn"fin'k'ir') (6) 

k^n^n' 

and treat it by perturbation theory. The first-order cor- 
rection to the eigenenergies Enk is given by 

SEl^k = J 'i'rj d'r'u:^{r)A{r,r')vnk{r') + h.c. (7) 

Inserting the eigenfunctions Mnfe(r) and w„fe(r) within 
the decoupling approximation, as given by eq. (2), and 
using the orthogonality if the Bloch wave- functions, 
!^^f'PnkiT)^n'k'{r) = dnn', One finds 

^^nk = '^nk'^nk ^ Afe_„j (r') (5„„i (5„2n + /l-C , 

= 0. (8) 

Eq. (3) allows us to interpret 2A„fc as the supercon- 
ducting gap. As a matter of principle, there is no guaran- 
tee that this gap reproduces the true gap of the supercon- 
ductor: Even with the exact exchange-correlation func- 
tional, the Kohn-Sham system need not reproduce the 



true spectrum of the fully interacting system. In semi- 
conductors and insulators, the fundamental gap is given 
by the Kohn-Sham gap plus the discontinuity of the xc 
potential with respect to the particle number. Standard 
functionals such as LDA and the GGAs are continuous 
with respect to the particle number and, therefore, can- 
not reproduce the discontinuity. In the superconducting 
case, the appearance of a discontinuity is not expected 
because wc are not working with fixed particle number 
in the first place. This fact alone does, of course, not 
prove that, for superconductors, the Kohn-Sham gap re- 
sulting from the exact xc functional would be identical 
with the true gap. This question remains the subject of 
future investigations. 



III. FUNCTIONALS 

A. fe-dependent functionals 

In this work we use the partially linearized versions 
for the functionals Znk and JCnk.n'k' proposed in I. The 
contributions stemming from the electron-phonon inter- 
action are obtained with the help of Kohn-Sham per- 
turbation theory. There are two terms: i) the first is 
non-diagonal 
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ph 

nk,n'k' 



Enk.n'k' 



tanh i^^^nk) tanh ^„'fe') a,, 

X [I{^nk,in'k' ,^\,q) - I{ink,-^n'k' ,^\,q)] , (9) 

where 5" q'" are the electron-phonon coupling con- 
stants and the function I is defined as 



(10) 



In the previous expression fp and Ufj are the Fermi-Dirac 
and Bose-Einstein distributions; ii) The second contribu- 
tion is diagonal in nk and reads 



^nk 



1 



\ nk.n'k' 



tanh y^ink j n'k' \,q 
[J{(,nk,in'k' ,^\,q) + J{^nk, -^n'k' ,^X,q)] 

where the function J is defined by 

j(e,^',n) = j(c,e',o)- j(e,^',-o). 

Finally we have 



(11) 

(12) 



(13) 
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On the other hand, the Coulomb interaction leads to the 
term 



K. 



with the definition 



TF 



■TF-ME _ TF 
nk.n'k' ^nk.n'k' ' 



(14) 



We then insert the unity 1 = /d^ — £,n'k') under the 
n'fc' summation in the right-hand side of Eq.(5), multiply 
the whole equation with the factor — Cnfe) ^'^d per- 
form the summation over all nk. Finally we divide the 
resulting equation by the density of states at the energy 
^, ^{0 = Enfc ^(^ - Cnfc)- This yields the gap equation 
in energy space, which now is only an one-dimensional 
integral equation: 



^*nk{rWnk{r')^n'k'{r)^n'k'ir'). (15) A(0 = -^(OA(0 



The electronic contribution is written in terms of the ma- 
trix elements (ME) of the screened Coulomb potential, 
as the use of the bare potential would be unrealistic. In 
this work we use a very simple model for the screening, 
namely the Thomas-Fermi (TF) model 



-kTF \r — r I 



\r — r' 



(16) 



with the Thomas- Fermi screening length, /ctfj given by 

= 47rAr(0) . (17) 

Finally, A^(0) denotes the total density of states at the 
Fermi level. Within this approach, the ME are calculated 
using the Bloch functions of the real material. In this 
basis, the ME read (using standard notation) 



TF _ \ ^ 

'^nk,n'k' / j 



(18) 



V ^ \k-k' + G\'^ + 

G 

X (n'fc'|e-'(''-'='+^)-'-|nfc) 
where V is the volume of the unit cell. 

B. Energy averaged functionals 



As discussed in the previous section, to obtain the gap 
function A„fe, we solve the gap equation (5) with the 
functionals given by Eqs. (9) and (11) for the electron- 
phonon interaction, and by Eq. (15) for the Coulomb 
repulsion. The inputs required for such calculation are 
the electron-phonon coupling constants 5" , and the 
normal-state Kohn-Sham eigcnencrgies i^„k and eigen- 
functions (fink- The coupling constants can be obtained 
from linear response^^ DFT calculations, while the eigen- 
states are the basic output of any standard DFT code. 

It is possible, however, to further simplify the solu- 
tion of Eq. (5). Very often, in the context of Eliashberg 
theory one neglects the gap anisotropy over the Fermi 
surface. We can apply the same approximation in our 
context (see Ref. 20 for further details). We assume that 
the pair potential is constant on iso-energy surfaces, i.e. 
it depends on the energies only: 



Anfe = A(^„fc) . 



(19) 



1 roc tanh (^E') 

-j d^'N{am,a — h^^io 



E' 



(20) 



The energy-averaged functions /C(^, ^') and Z{(^) are de- 
fined as: 



1 



5{£,-£,nk)^{i' -£,n'k')Knk,n'k' 



nk,n' k' 



- ^nk)Znk ■ 



(21) 

(22) 



nk 



We should mention that in performing the average over 
iso-energetic surfaces we assumed an s-wave pairing field. 
It is straightforward to devise similar averaging proce- 
dures for pairing fields of different symmetry. 

The phononic contributions to the averaged function- 
als read 



/cp^(^,r) 



1 



tanh (fe) tanh(fe) ^(0) 



dfi a^Fin) 



(23) 



and 



1 



DO /> 



tanh (1^) 

X [j(^,^',o) + j(^,-^',n)] 



(24) 



with / and J given by Eqs. (10) and (12). Finally, the 
Eliashberg spectral function is the electron-phonon cou- 
pling constant averaged on the Fermi surface 



1 



iV(0) 



nk,n'k' 



nk,n'k' X,q 

xSi^nk)5{Cn'k')S{n-nx,g). (25) 



Note that in Eq. (25) (and, consequently, in Eq. (23)) 
we replaced the density of states N{$,) by its value at 
the Fermi energy A'^(O). This procediue is well justified, 
because it only requires that for each given band the cou- 
plings do not change much on an energy scale of the order 
of the debye frequency, i.e., meV. The energy dependence 
of the density of states, on the other hand, is kept in Eqs. 
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(20) and (22), as N{(_) can vary even on this small energy 
scale (e.g., in transition metals where the Fermi energy 
is at the edge of large peaks in the density of states). 

In order to evaluate the Coulomb terms, we further ap- 
proximate the Kohn-Sham eigenvalues by a free-electron 
(FE) parabolic dispersion Sk — fc^/2. This approxi- 
mation is well justified for simple metals, but it is ex- 
pected to fail for more complicated systems. Within this 
approximation the energy-average of the Thomas-Fermi 
screened Coulomb interaction, given in Eq. (18), reads: 



TF-FE 



(26) 



with k = ■y/2{^ - n) and k' = ■s/2{£,' - fj,). 



C. Hybrid functionals 

It is possible to obtain a hybrid approach where the 
averaged functionals are used in the fc-dependent gap- 
equation (5). The averaged phononic terms (23) and (24) 
can be used in (5), just by replacing the energies ^ and 
^' by the energy eigenvalues of the real material ^nk and 



Zfk = 



(27) 
(28) 



For the electronic terms we use an approach that goes 
along the lines of Sham and Kohn^^ (SK), as detailed 
in I. The basic idea is again to replace the free-electron 
bands by the real bands of the material, and furthermore 
to adjust the chemical potentials of the two systems. The 
functional reads 



j^-TF-SK _ ^ 
'^nk,n'k' o /Z^ ^ 



X log I 



Vnk + Vn'k' + 2^ 
J]nh H" '^n'h' ~~ ^/ Vnk'Hn' k' + ^Tf/^> 



(29) 



In this expression we defined rjnk = S,nk + ^^F' where fcp 
is the Fermi wave- vector of a homogeneous electron gas 
with (constant) density equal to the average density of 
the material. 



IV. COMPUTATIONAL DETAILS 

Obtaining superconducting properties through the so- 
lution of the gap equation can be viewed as post- 
processing results of standard electronic structure cal- 
culations. In this work we used electronic band struc- 
tures obtained from full-potential linearized augmented 
plane wave^^ calculations; the same method is also used 
to compute the Coulomb potential matrix elements (18), 
along the lines described in Ref. 23. On the other hand. 



the Eliashbcrg function can be obtained from linear re- 
sponse calculations^^. (Note that the Eliashberg func- 
tions can also be extracted from experiment. While the 
use of such experimental a^F^s renders the theory semi- 
phenomenological it often gives useful insights.) 

We developed two completely independent codes: the 
first solves the averaged gap equation (20) in energy 
space; the second solves the fc-resolved gap equation (5). 
The two schemes were compared whenever possible, and 
give similar results for the simple metals studied here (see 
the discussion at the end of this section). The solution 
of the averaged gap equation is a quite simple numerical 
task. The energy is discretizcd in a logarithmic mesh in 
order to increase the accuracy close to the Fermi level 
and the resulting equation is iterated until convergence. 
This code was used within the TF-FE scheme. (For a 
summary of the difFca'ent schemes, please refer to Ta- 
ble I.) On the other hand, to use the TF-ME and the 
TF-SK functionals we need to solve the fe-resolved equa- 
tion. This is a difficult numerical task mainly because: i) 
we need to describe the energy bands over a large energy 
window and, at the same time, ii) we require a very good 
resolution in a small energy window around the Fermi en- 
ergy. Requirement i) comes from the large energy scales 
involved in the Thomas-Fermi screened Coulomb poten- 
tial (typically of the order of the Fermi energy). This 
is not an artifact of the static nature of the Thomas- 
Fermi screening. In fact, even a more realistic dynami- 
cally screened potential would need a very large energy 
range. This point will be further discussed in the fol- 
lowing section. To illustrate requirement ii) we plot, in 
Fig. (1), the kernel of Eq. (23), 

W{^, r, O) = J($, f , fi) - J(^, -e, fi) . (30) 

The left panel of Fig. 1 depicts the variation of W with 
energy. This plot was obtained at very low tempera- 
ture T = 0.01 K, by fixing ^' very close to the Fermi 
surface ^ « 0. (We did not use the values T = and 
^ = to avoid numerical problems.) Considering that 
the kernel decreases to about to about 20% of its value 
within roughly 10 meV, it is clear that we need to sample 
very carefully this energy range. This is particularly im- 
portant for materials, like niobium and tantalum, where 
the Fermi energy is at the edge of a peak in the den- 
sity of states. In the right panel we plot the dependence 
of W on the frequency fl, for ^ and ^' very close to the 
Fermi surface. The kernel W acts as a weight for a^F{ft) 
(see, e.g., Eq. (23)), and enhances the contribution of the 
lower-frequency phonons. 

In order to fulfill the requirements i) and ii) we devel- 
oped the following numerical framework: Starting from 
ab-initio bands calculated over a few hundred fc-points 
in the irreducible wedge of the Brillouin zone (BZ), we 
compute a good spline fit of S^nk over a Fourier series, 
according to the scheme of Koelling and Wood^^. Using 
this fit, we then obtain the energies ^ over a very large set 
of random fe-points, distributed according to a Metropo- 
lis algorithm. This algorithm was devised to accumulate 
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^' [eV] Q. [THz] 

FIG. 1: Universal kernel W{(,,^' ,0.) as a function of ^' and Q. Left panel: universal kernel W{(, ~ 0,^',n = ITHz) as a 
function of 5'. Right panel: = ^' « 0, fl) as a function of fl. Both panels were obtained for T = 0.01 K. 



a largo number of fe-points in the first few meV's around 
the Fermi surface. A good convergence can be reached 
by using about 15000-20000 independent points for each 
band crossing the Fermi surface, while a reasonable de- 
scription of the remaining bands is obtained with about 
1000 independent points per band. The TF-ME scheme 
is the slowest to converge, as it involves the integration 
of a fimction which is peaked around |fc' — fc| (the 
integrand would diverge if we used the bare Coulomb 
potential instead of the Thomas-Fermi screened interac- 
tion). 

The function A„fc only needs to be defined in the irre- 
ducible wedge of the Brillouin zone. On the other hand, 
when using the TF-ME scheme, the k' summation has to 
be performed over the whole zone. This can be seen from 
the expression of the Coidonib matrix elements, Eq. (18): 
By setting fe ^ the symmetry of the k' summation 
is broken (only the operations of the little group of k 
can be retained). This is a well-known situation also 
in the framework of electronic self-energy calculations. 
The case of the hybrid functional is simpler. As these 
functionals depend on k and k' though ^nk and (,n'k', 
and as the eigenvalues are totally symmetric with re- 
spect to the crystal point group, the summation over k' 
can be limited to the irreducible wedge of the Brillouin 
zone. This formalism can be easily generalized to the case 
of non .s-wave pairing, by imposing that the gap trans- 
forms according to a given irreducible representation of 
the crystal point group. 

As an initial guess for the gap function A„fc we use a 
step function. By using a Broyden scheme^^ to mix Ay^^, 
we obtain convergence in a mere 10-15 self-consistent it- 
erations. The converged result at a given temperature 
can be used as a starting point for the next tempera- 
ture. Clearly, this procedure reduces the total number of 
iterations required. 

The matrix elements of the screened Coulomb poten- 
tial are obtained with the full-potential linearized aug- 
mented plane wave method, as explained in Ref. 23. 
They are calculated for all bands within 15-30 eV from 



the Fermi surface (typically 14 valence and conduction 
bands). It is unfeasible to obtain the matrix elements for 
the whole set of random fe-points used in the solution of 
the gap equation. Thus, we calculate the matrix elements 
for a set of fc and k' belonging to a regular mesh. The 
values at the random points k and fe' are then obtained 
by mapping fe and fe' to the corresponding hypercube of 
the regular mesh. Test calciilations have shown that al- 
ready a 6 X 6 X 6 mesh gives acceptably small residual 
errors (1-2 % of the gap at Ep). A finer mesh is however 
needed to quantitatively estimate the spread of the gap 
for each given energy, which is material dependent and 
mostly due to the physical fe and fe' dependence of the 
matrix elements. 

In Fig. 2 we show a detailed analysis of the conver- 
gence of our method for Nb, using the TF-SK approach. 
In this figure we plot the relative differences between the 
k-resolved approach and the results of the energy-only 
scheme, which can be considered to be numerically con- 
vergent (because of its use of a logarithmic energy mesh), 
but itself depending on the fine details of the input den- 
sity of states. This is particularly true in the case of Nb, 
where Ep sits in a rapidly varying part of the DOS. The 
results are shown as a function of the number of inde- 
pendent k-points used for each band crossing the Fermi 
level, and circles and squares represent the results ob- 
tained using two different distribution functions used to 
sample the Brillouin zone (as detailed in the figure cap- 
tion). We can see the capability of our procedure to 
converge towards the (computationally completely inde- 
pendent) energy-only result, with a very small numerical 
error. The two rather different samplings lead to almost 
identical averages, with a correct, gaussian-like, distribu- 
tion of results around the average. We notice that bands 
away from the Ep are sampled uniformly in k-space, as 
it should since at these cnc^rgies only the Coulomb term 
remains, with a weak energy dependence. 

For simplicity, we used in all our calculations the av- 
eraged or hybrid phononic functionals. In this case, 
the electron-phonon interaction enters the calculation 
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FIG. 2: (Color online) Relative difference between the k- 
resolved and the energy-integrated results for Nb, as a func- 
tion of the number of independent fe-points used for each band 
crossing Ef- The circles and squares are obtained, respec- 
tively, with probability distributions equal to max{x{£,, T — 
0), 0.05} and to a stepwise function having widths of 2 and 40 
mRy, as shown in the insert. Thin lines represent the arith- 
metic averages of the k-resolved data. 

through the EHashberg function a^F{il). This cor- 
responds to neglecting the anisotropy of the electron- 
phonon coupling, which should be a good approximation 
for the particular systems studied in the following. 

We quantify the precision of our results to be around 
5%. This value includes the uncertainty associated with 
the Metropolis procedure. We emphasize that this is 
a non-trivial numerical achievement, as superconducting 
properties depend exponentially on the electron-phonon 
coupling and on the Coulomb interaction. Particularly 
difficult is to obtain numerically stable results for the 
small gap materials, where superconductivity stems from 
an almost complete cancellation of large and opposite 
terms. 



V. RESULTS AND DISCUSSION 

In the following we will present results obtained for the 
simple metals Al, Mo, Ta, Nb, and Pb. Note that this 
group of materials includes both weak couphng (Al, Mo, 
and Ta) and strong coupling (Nb and Pb) superconduc- 
tors. We solved the gap equation within three different 
approaches, that we label TF-FE, TF-SK, and TF-ME 
(see Table I). 

In Fig. 3 we show the Eliashberg functions used in 
our calculations. All these curves were obtained by 
Savrasov"^^ using linear response theory. The five ma- 
terials cover a broad frequency range, from 1 to 10 THz. 
The a^F{^) function for Pb is located at much lower fre- 
quencies than for the other materials, due to the heavy 
nuclear mass of Pb. For this reason, Pb has the largest 
total electron-phonon coupling constant A of the materi- 
als studied (A = 2 JdQ a^F{fl)/fl). Nevertheless, Nb be- 



FIG. 3: (Color online) The Ehashberg function a'^F{n) for 
the simple metals used in our calculations. The curves are 
taken from Ref. 31. 



comes superconducting at higher temperatures than Pb. 
This can be easily understood on the basis of McMillan's 
formula:'^^ Besides the well-know exponential dependence 
on the electron-phonon coupling constant, the transition 
temperature increases linearly with the average phonon 
frequency. Lowering the phonon frequencies has therefore 
both a positive and a negative effect on Tc. An accurate 
description of the superconducting state must take into 
account both these effects. 

Before examining our results, we depict in Fig. 4 the 
behavior of the phononic functionals, —N{0)ICp^ 
and ZP^(^). Both terms turn out to be very peaked at 
the Fermi energy, and are non-zero only in a small region 
around it (of the order of the Debye energy). Further- 
more, the values of the functionals at the Fermi energy 
read 7V(0)/Cp'^(0, 0) = -A, and Zp'^(O) = A. The K.^^ 
term is negative, reflecting the electron-phonon attrac- 
tion responsible for the superconducting state. On the 
other hand, the term Zp'^ is positive and tends to de- 
crease the superconducting gap and therefore Tc. 

In Fig. 5 we show the superconducting gap A calcu- 
lated using different approaches for Pb and Nb. For Pb 
we also show curves calculated at different temperatures. 
It is important to note that for T > Tc the self-consistent 
calculation correctly converges to zero even if the starting 
gap function is different from zero. The curves labeled 
TF-ME show the detailed sampling of the whole energy 
range given by our fe-point mesh, down to a few fieV. 
Within the TF-ME approach, the spread in the values 
of the gap function of Pb at energies near Ep is quite 
small (around 2%), which can be understood as a conse- 
quence of a rather isotropic Fermi surface. This spread is 
larger for the higher energy states because of the larger 
spread of orbital character in the electronic states. The 
gap function for both materials exhibits a very similar 
shape, with a node at about 0.4 meV, and a negative 
tail extending to high energy. This shape is basically in- 
dependent of the temperature, and is a general feature 
found in all our calculations. Due to the presence of 
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Method Gap-equation 



Phononic terms 



Coulomb term 



TF-ME 



TF-SK 



TF-FE 



fe-dependent 

Eq. (5) 
fe-dependent 

Eq. (5) 

Energy averaged 
Eq. (20) 



Hybrid: averaged with real bands 
Eqs. (23) and (24) 

Hybrid: averaged with real bands 
Eqs. (23) and (24) 

Averaged 
Eqs. (23) and (24) 



Matrix elements 
Eq. (18) 

Hybrid: averaged with real bands 
Eq. (29) 

Averaged 
Eq. (26) 



TABLE I: Summary of the different methods used in our calculations. 




FIG. 4: (Color online) Diagonal of ~N{0)1C^^{^, (left panel) and Z^^^i^) at T = 0.01 K for the simple metals studied. 



Coulomb repulsion, this shape is a necessary condition 
to obtain a superconducting solution. It is well known 
that the structure of the gap equation is such that the 
repulsive Coulomb interaction between two Cooper pairs 
(at {k Ij—k 1} and {k' 1,—k' |}) gives a constructive 
contribution to the gap if the values of A^k and A„/fc/ 
have opposite signs. This condition is realized when, e.g., 
^nfe is small and Cn'k' is large. Similar arguments are the 
basis of the classical calculations of Ref. 33, and result 
in the definition of the renormalized Coulomb parameter 
fi*. 

The three different schemes (TF-FE, TF-SK, and TF- 
ME) agree well among themselves for Pb, while for Nb 
the TF-SK curve differs by 10% from the TF-SK and TF- 
ME results. This difference can be explained by looking 
at the band structure of the materials. In Pb, the va- 
lence and conduction bands are basically due to s and p 
orbitals, for which the averaged and hybrid schemes work 
quite well. The same result is found for Al, with an even 
larger similarity between TF-FE, TF-SK and TF-ME re- 
sults. All this is easily understood from the fact that the 
three schemes TF-FE, TF-SK, and TF-ME become iden- 
tical in the limit of the imiform gas. Hence one would 
expect them to yield similar results for the delocalized s- 
p orbitals. On the other hand, the bands of Nb exhibit a 
strong d-character, which leads to the difference between 
the methods. A similar behavior is found for Ta. This 
trend is also confirmed by the value of the gap function 
for the low lying semi-core d states of Pb (not shown in 
the figure) which is very different in the TF-SK, and TF- 



ME calculations (it is much smaller for TF-ME) . This is a 
clear consequence of the localized nature of these states, 
implying a very large repulsive self-term for those bands. 
In principle, the TF-ME approach should be the most 
precise. 

Note that our approximate functionals use the 
Thomas-Fermi model for screening. We have compared 
the matrix elements of the Thomas-Fermi-screened po- 
tential with those coming from a complete (static) RPA 
and found very close agreement between the two ap- 
proaches. It is very important in this context to define 
the Thomas- Fermi wave vector kxF in terms of the den- 
sity of states at the Fermi level, N{0). As the latter is 
obtained from a full-scale Kohn-Sham calculation, kxF 
is treated here implicitly as a rather complicated density 
functional. 

Furthermore, we neglect the off-diagonal elements of 
the dielectric matrix. This is not necessarily a good ap- 
proximation for transition metal compounds. However, 
it was shown in the literature^^'^^ that the changes pro- 
duced by these local field effects are largely cancelled 
by the additional inclusion of exchange-correlation effects 
(generated by the exchange-correlation kernel^* of time- 
dependent density functional theory^^). We believe that 
our good results, obtained by neglecting both local-field 
and exchange-correlation effects, can at least partly be 
explained by this cancellation. Further help comes from 
the fact that superconductivity describes correlations on 
the scale of the coherence length, which is usually much 
larger than the dimensions of the unit cell. On this scale, 
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FIG. 5: (Color online) The function A(^„fc,T) for lead (left panel) and niobium (right panel). 



Tc [K] 





TF-ME 


TF-SK 


TF-FE 


exp 


A 


Mo 




0.33 


0.54 


0.92 


0.42 


Al 


0.90 


0.90 


1.0 


1.18 


0.44 


Ta 


3.7 


2.7 


4.8 


4.48 


0.84 


Pb 


6.9 


7.2 


6.8 


7.2 


1.62 


Nb 


9.5 


8.4 


9.4 


9.3 


1.18 






Ao 


[meV] 








TF-ME 


TF-SK 


TF-FE 


exp 


A 


Mo 




0.049 


0.099 




0.42 


Al 


0.14 


0.15 


0.15 


0.179 


0.44 


Ta 


0.63 


0.53 


0.76 


0.694 


0.84 


Pb 


1.34 


1.40 


1.31 


1.33 


1.62 


Nb 


1.74 


1.54 


1.79 


1.55 


1.18 



TABLE II: The critical temperature (upper panel) and the 
superconducting gap at Fermi level and T — 0.01 K (lower 
panel), compared with experiment'^". We show also the total 
electron-phonon coupling constant A''^. 



local field effects are certainly less important. 

The gap function at the Fermi energy is depicted in 
Fig. 6 as a function of temperature for Pb and Nb. The 
curves show a BCS-like square root behavior close to Tc, 
which we would expect for these simple superconductors. 
For Pb, both Tc and Aq agree quite well among them- 
selves and with the experimental results. For Nb, on 
the other hand, the TF-SK gap is roughly 15% smaller 
than the TF-ME curve. Curiously, the TF-SK gap at 
zero temperature is very close to the experimental value, 
while the TF-ME yields a much better transition tem- 
perature. The difference between the methods can be 
once more explained by the d-character of the bands in 
Nb close to the Fermi energy. 

Our results for the superconducting transition temper- 
atures and gap functions are summarized in Table II and 
in Fig. 7 for all materials studied. In the same table 
we also show the experimental results, and the values of 
the electron-phonon coupling constant A. Mo is a weak 
coupling superconductor with a very small gap and very 



low transition temperature. In this case, the different 
theoretical approaches lead to quite different results, and 
the overall agreement with experiment is not very good. 
Surprisingly, it is the TF-FE approach which gives the 
better results. However, many tests showed that the Mo 
results are very sensitive to the details of the density 
of states underlying our calculations. This is quite nor- 
mal for a material with such a small gap resulting from 
a very fine balance between large terms. In Al the elec- 
tronic states have a strong free-electron character and the 
density of states follows closely the free-electron DOS. It 
is not surprising, therefore, that all three methods yield 
very similar results. These results are also in satisfactory 
agreement with experiment. Ta shows a behavior simi- 
lar to Nb, but the agreement with experiment is poorer. 
Moreover, the larger spread of values for Ta relative to 
Nb can be explained by the smaller values of A and Tc- 
As in the case of Mo, the TF-FE approach yields the 
best results. This surprising result can be understood to 
some extent: The TF-FE approach is fully consistent, in 
the sense that all the quantities entering the method are 
averaged in similar ways; In the TF-SK approach, on the 
other hand, certain quantities (but not all) are calculated 
at the average density. We emphasize, however, that the 
agreement of our results with experiment, without mak- 
ing use of any adjustable parameters, is unprecedented 
in the field. We note in passing that for Cu we could 
not find a non-trivial solution of the gap equation at any 
temperature, which also in in agreement with the exper- 
imentally observed absence of superconductivity in Cu. 

It is instructive to visualize how the different terms en- 
tering the gap equation combine to yield the values listed 
in Tab. II: The phononic term K.^^ is negative, and is the 
responsible for the superconducting state in the simple 
metals we studied; Z^^ gives a repulsive contribution, 
whose relative importance increases for the strong cou- 
pling materials; finally, the electronic Coulomb repulsion 
leads to a large cancellation between constructive and de- 
structive interference effects. To better understand these 
effects we plot, in Fig. 8, the gap function at zero tem- 
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FIG. 6: (Color online) The gap function at the Fermi surface Aq for Pb and Nb as a function of temperature. 




FIG. 7: (Color online) The critical temperature and, superconducting gap at Fermi level and T = 0.01 K, compared with 
experiment. The numerical values can be found in Table II. 



perature of Pb obtained through the solution of the gap 
equation (20) including only /Cp^, the two phononic con- 
tributions, or all three phononic and Coulomb contribu- 
tions. Clearly, the value of the gap is the result of a 
subtle interplay of opposite contributions, each one con- 
siderably larger than the gap itself. 

Another interesting point regards the convergence of 
the gap function with the energy cutoff. This conver- 
gence can be observed in Fig. 9, where we plot the gap for 
Pb and Al calculated within the TF-FE approach. The 
different materials behave differently, with a faster con- 
vergence for the material with stronger elect ron-phonon 
coupling. However, even for Pb, an energy cutoff of at 
least 10 eV was necessary to achieve convergence. It is 
a key feature of our approach that the matrix elements 
of the screened Coulomb interaction are used in the ker- 
nel of the gap equation up to very high energies. It is 
this feature that allows for the description of non-trivial, 
material-specific effects. Traditionally, the use of this 
large energy window is avoided by rescaling the Fermi- 



surface average, ^, of these matrix elements. This rescal- 
ing leads to the the Morel- Anderson pseudopotential ^* . 
While ^* is an ingenious concept to capture the essential 
physics of the Coulomb repulsion in simple superconduc- 
tors, it is not sufficient to treat more complex materi- 
als such as, e.g., MgB2. The detailed analysis carried 
out in Ref. 5 shows how the actual value of matrix el- 
ements of the screened Coulomb interaction, computed 
w.r.t. Bloch states of different orbital character (ct and 
tt), leads to non-trivial effects, determined by the specific 
character of the orbitals. This indicates that, contrary to 
common wisdom, superconducting properties are not ex- 
clusively determined by a small region around the Fermi 
level. Only the contributions from states higher than 
20-30 eV are essentially negligible, due to the asymptotic 
decay of the Coulomb term (26). 



After solving the gap equation, it is straightforward 
to calculate thermodynamic functions, such as the elec- 
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FIG. 8: (Color online) Gap function versus energy at T = 
0.01 K for Pb, including the different contributions to the gap 
equation. 

Theory Experiment 



Pb 2.93 

Nb 2.87 

Ta 2.64 

Al 2.46 



3.57-3.71 
2.8-3.07 

2.63 

2.43 



TABLE III: 



Normalized electronic specific heat 
as computed from the TF-ME approach. 



tronic entropy of the Kohn-Sham system 

= -2fcB {/M^nfe) log [fl3{Enk)] 
nk 

+ [1 - fp{Enk)] log [1 - l0{E^k)] } , (31) 

where is the Bohzmann constant. In Fig. 10 we de- 
pict the difference of entropy between the superconduct- 
ing and the normal states A^* = S's — for Pb and 
Nb. This quantity has the expected temperature depen- 
dence, going smoothly to zero at Tc. This indicates the 
stability of our calculations even close to the transition 
temperature where the gap becomes very small. In the 
same figure we also depict the normalized specific heat 
Cf'^(r)/CN(rc). The specific heat is obtained by eval- 
uating numerically the temperature derivative of the en- 
tropy. Despite the numerical uncertainty associated to 
the calculation of the derivative, the curves of the spe- 
cific heat are quite stable. The discontinuities of the spe- 
cific heat at obtained within the TF-ME approach are 
shown in Table III. Our results are in quite good agree- 
ment with experiment: While for Al and Ta we confirm 
the BCS value found in experiments, we reproduce the 
strong coupling value of Nb. For Pb the agreement with 
experiment is worse, but still within acceptable margins. 

It is clear that any theory that aims at describing the 
superconducting state has to include retardation effects. 
While the main equation of our density functional formal- 



ism, the gap equation (5), has the form of a static equa- 
tion, retardation enters through the electron-phonon part 
of the exchange-correlation potential. To further demon- 
strate that our theory takes retardation effects properly 
into account, we calculated the isotope effect coefficient 
a. In fact, if retardation effects were not included, a 
would be equal to the BCS value a — 0.5. For a mono- 
atomic solid, the Eliashberg function a^F{Q) for mate- 
rials with different isotopic masses can be obtained by 
rescaling the phonon frequencies with the square root 
of the nuclear mass M . Note that only the frequen- 
cies are rescaled, while the total electron-phonon cou- 
pling constant A remains unchanged. Within the TF-FE 
approach, we computed the isotope effect of the gap func- 
tion. We obtained the values 0.37 and 0.47 for Mo and Pb 
respectively, which can be compared to the correspond- 
ing experimental values, a = 0.33 and a — 0.47. This 
agreement with experiment, in the presence a significant 
deviation from the BCS value, proves, in our opinion, 
that retardation effects are properly taken into account 
in our calculations. 

To conclude the presentation of our results we show, 
in Fig. 11 the order parameter represented on the basis 
of Bloch orbitals. The lower panel of the figure depicts 
the order parameter for Al, Pb and Nb as a function 
of the energy of the corresponding Bloch state. As ex- 
pected, Xnk is very localized in a small energy window 
around the Fermi level. The spread of the order param- 
eter in fc-space is related to the inverse of the coherence 
length of the superconductor. This can be visualized by 
rescaling the relative energy ^ by S^o/hv-p, where is 
the experimental coherence length of the material and 
hvp = Vfc^nfc is the Fermi velocity. The resulting curves 
depict the spread in fc-space of the wave-packets describ- 
ing the Cooper pairs in units of the experimental coher- 
ence length. As expected, the curves for Al, Pb, and Nb 
are quite similar to each other, and have a width com- 
parable to unity. We can therefore conclude that not 
only the maximum gap value but also its overall energy 
dependence are described accurately in our approach. 



VI. CONCLUSIONS 

In this work we present the first full-scale applica- 
tion of the ab-initio theory for superconductivity, which 
was developed in the preceeding paper (I). Superconduct- 
ing properties of simple conventional superconductors are 
computed without any experimental input. In this way, 
we were able to test the theory developed in I and to 
assess the quality of the functionals proposed. It turns 
out that the different proposed functionals lead to re- 
sults which are in good agreement with each other for 
the simple metals studied. The agreement is better when 
the normal-state Bloch functions are not too strongly lo- 
calized (as, e.g., the sp-orbitals in metals in contrast to 
the more localized d states of the transition elements). 
The most important result is that the calculated transi- 
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FIG. 9: (Color online) Convergence of the calculated T — 0.01 K gap as a function of the energy cutoff, for Pb and Al. 
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FIG. 10: (Color online) Thermodynamic functions for Pb (left panels) and Nb (right panels). Upper panels: ratio between the 
electronic specific heat ratio in the normal (N) or superconducting (S) states and the specific heat at Tc (Cf '^(T)/C^(Tc)). 
Lower panels: difference of entropy between the superconducting and the normal states AS* = 5*3 — Sn- 




^-2 -1.5 -1 -0.5 0.5 1 1.5 2 



0.5 
W 0.4 



— 1 — 




-1 — ' — 1 


1 — ■ — 1 — ' — 1 — ■ — 1 — 

i 






i 
1 

i 
i 


i : 
t 






i 













_I , ■ ■ • . -^-^ - - — ■ . 1_ 

-0.06 -0.04 -0.02 0.02 0.04 0.06 
^JeV] 



FIG. 11: (Color online) Order parameter for Al, Pb, and Nb 
represented on the basis of Bloch orbitals as a function of 
momentum (upper panel) and energy (lower panel). The mo- 
mentum was rescaled with the experimental coherence lengths 
^0 = 1600, 90, and 40 nm for Al, Pb and Nb respectively. 



tion temperatures and superconducting gaps are also in 
good agreement with experimental values. The largest 
deviations from the experimental results are found for 
the elements in the weak coupling limit with Mo being 
the most pronounced example. We emphasize, however, 
that the agreement of our results with experiment, with- 
out making use of any adjustable parameter, is unprece- 
dented in the field. Furthermore, we also obtained other 
quantities such as the entropy and the specific heat. Our 
approach reproduces the correct values for the disconti- 
nuity of the specific heat at even in the strong cou- 
pling regime. Finally, we calculated the isotope effect 
for Mo and Pb, achieving again rather good agreement 
with experiment. These results clearly show that retar- 
dation effects are correctly described by the theory. Our 
calculations demonstrate that, as far as conventional su- 
perconductivity is concerned, the ab-initio prediction of 
superconducting properties is feasible. 
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